      double precision Function EOS(x,key)
      use accur
      Implicit none
      double precision::x,rho,rho_,rho_c,p,p_c,gamma,CK,H,F
      integer::key,i,IMax,IMax0,io
      Parameter(IMax0=500)
      Dimension rho(0:IMax0),p(0:IMax0),gamma(IMax0+1),CK(IMax0+1),
     *	H(0:IMax0),F(0:IMax0)
      Logical LogInit
      Data LogInit/.true./
      if(LogInit)then
C      open(3,file='EOS_s15.dat')
!       open(3,file='Polytrope_n1.dat')
      open(3,file='Potekhin19.dat')
      rho(0)=zero
      p(0)=zero
      IMax=0
      do i=1,IMax0
!      read(3,*,err=100)rho_c,p_c
      read(3,*,iostat=io)rho_c,p_c
      if(io<0)exit;
      if(rho_c.gt.zero)then
      IMax=IMax+1
      rho(IMax)=rho_c
      p(IMax)=p_c
      endif
      enddo
100   close(3)
      do i=2,IMax
      gamma(i)=(log(p(i))-log(p(i-1)))/(log(rho(i))-log(rho(i-1)))
      CK(i)=exp((log(p(i-1))*log(rho(i))-log(p(i))*log(rho(i-1)))/
     *	(log(rho(i))-log(rho(i-1))))
      enddo
      gamma(1)=gamma(2)
      gamma(IMax+1)=gamma(IMax)
      CK(1)=CK(2)
      CK(IMax+1)=CK(IMax)
      H(0)=zero
      H(1)=CK(1)*gamma(1)*(rho(1)**(gamma(1)-one))/(gamma(1)-one)
      do i=2,IMax
      H(i)=H(i-1)+CK(i)*gamma(i)*(rho(i)**(gamma(i)-one)-rho(i-1)**
     *	(gamma(i)-one))/(gamma(i)-one)
      enddo
      F(0)=zero
      F(1)=CK(1)*(rho(1)**(gamma(1)-one))/(gamma(1)-one)
      do i=2,IMax
      F(i)=F(i-1)+CK(i)*(rho(i)**(gamma(i)-one)-rho(i-1)**
     * 	(gamma(i)-one))/(gamma(i)-one)
      enddo
      LogInit=.false.

      endif
      if(key.eq.1)then
      if(x.gt.zero)then
      i=IMax+1
      do while(x.lt.H(i-1))
      i=i-1
      enddo
      if(i.eq.1)then
      EOS=(x*(gamma(i)-one)/(gamma(i)*CK(i)))**(one/(gamma(i)-one))
      else
      EOS=((x-H(i-1))*(gamma(i)-one)/(gamma(i)*CK(i))+rho(i-1)**
     *	(gamma(i)-one))**(one/(gamma(i)-one))
      endif
      else
      EOS=zero
      endif
      endif
      if(key.eq.2)then
      if(x.gt.zero)then
      i=IMax+1
      do while(x.lt.H(i-1))
      i=i-1
      enddo
      if(i.eq.1)then
      EOS=(x*(gamma(i)-one)/(gamma(i)*CK(i)))**(one/(gamma(i)-one)-one)
     *	/(CK(i)*gamma(i))
      else
      EOS=((x-H(i-1))*(gamma(i)-one)/(gamma(i)*CK(i))+rho(i-1)**
     *	(gamma(i)-one))**(one/(gamma(i)-one)-one)/(CK(i)*gamma(i))
      endif
      else
      EOS=zero
      endif
      endif
      if(key.eq.3)then
      if(x.gt.zero)then
      i=IMax+1
      do while(x.lt.rho(i-1))
      i=i-1
      enddo
      if(i.eq.1)then
      EOS=H(i-1)+CK(i)*gamma(i)*(x**(gamma(i)-one))/(gamma(i)-one)
      else
      EOS=H(i-1)+CK(i)*gamma(i)*(x**(gamma(i)-one)-rho(i-1)**
     *	(gamma(i)-one))/(gamma(i)-one)
      endif
      else
      EOS=zero
      endif
      endif
      if(key.eq.4)then
      if(x.gt.zero)then
      i=IMax+1
      do while(x.lt.H(i-1))
      i=i-1
      enddo
      if(i.eq.1)then
      rho_=(x*(gamma(i)-one)/(gamma(i)*CK(i)))**(one/(gamma(i)-one))
      else
      rho_=((x-H(i-1))*(gamma(i)-one)/(gamma(i)*CK(i))+rho(i-1)**
     *	(gamma(i)-one))**(one/(gamma(i)-one))
      endif
      EOS=CK(i)*rho_**gamma(i)
      else
      EOS=zero
      endif
      endif
      if(key.eq.5)then
      if(x.gt.zero)then
      i=IMax+1
      do while(x.lt.rho(i-1))
      i=i-1
      enddo
      if(i.eq.1)then
      EOS=F(i-1)+CK(i)*(x**(gamma(i)-one))/(gamma(i)-one)
      else
      EOS=F(i-1)+CK(i)*(x**(gamma(i)-one)-rho(i-1)**(gamma(i)-one))
     *	/(gamma(i)-one)
      endif
      EOS=EOS*x
      else
      EOS=zero
      endif
      endif
      return
      end Function EOS
